/* ======Trend Inflation and the Nature of Structural Breaks in the New Keynesian Phillips Curve ======

This program estimates the unobserved components model for the NKPC with stochastic trend inflation
with two structural breaks. The variability of shocks to trend inflation are constrained to be constant
across regimes.

Written by: Pym Manopimoke

Model:

Pi_t = Pi_t^* + k_{S_t} Sum_{j=0}^{inf} beta^j E_{t-1}(x_{t+j}) + z_t
Pi_t^* = Pi_{t-1}^* + v_t, v_t ~iid.N(0,sig^2_v)
x_t = phi_1 x_{t-1} + phi_2 x_{t-2} + u_t, u_t|D_t~iid.N(0,sig^2_{u,D_t})
z_t = psi_{S_t} z_{t-1} + epsilon_t, epsilon_t|S_t ~iid.N(0,sig^2_{epsilon,S_t})

u_t and epsilon_t are correlated with correlation coefficient rho.

S_t = {0,1}
Prob[ S(t)=0 | S(t-1)=0] = p00
Prob[ S(t)=1 | S(t-1)=1] = 1

sig^2_{u,D_t} = sig^2_{u,0}(1-D_t) + sig^2_{u,1}D_t
D_t = {0 if t<= 1984Q3, 1 otherwise}

========================================================================================================*/

new;
library optmum,PGRAPH;
pqgwin auto;
cls;
 
load data[250,11]=data.prn;  
data=data[2:250,.];
t=rows(data);
@Set the sample to 1952:1-2007:3@
startsmpl=21;
endsmpl =243; 

qtr= data[startsmpl:endsmpl,1];  
/*Column 2 for GDP chained-weighted inflation, 7 for CPI inflation, 8 for PCE inflation*/  
inf=data[startsmpl:endsmpl,8];  @PCE inflation@
inff1=data[startsmpl+1:endsmpl+1,8];  @one-quarter ahead inflation@

xx=data[startsmpl:endsmpl,5];  @CBO output gap@
t=rows(inf~xx);

@Define dummy variable series D_t to account for the Great Moderation@
GM = 130; @1984:2@
dmat = ones(GM,1)|zeros(t-GM,1); 

@================= Initialize Global Variables============@

tvp_no=4; 
start= 11;
t1=75;t2=125;

prmtr_in=   4.3360830 |
            1.0297822 |
            0.45991641 |
            0.34008988 |
            0.73952899 |
            1.5544761 |
            -0.083434780 |
            1.4732290 |
            -4.3488046 |
            -0.077438017 |
            9.2330093 |
            0.018851190|
            0 ;
{xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN);

PRM_FNL=TRANS(xout);             @ Estimated  coefficients, constrained@
xout;
hout0=hessp(&lik_fcn,xout);
hout0;
hout=inv(hout0);
grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl)); @Standard errors of the estimated coefficients@
 
{B_TT,PR_TT,PR_TL,stdprederr,xfor,confint} =FILTER(XOUT);

_pdate=0;
_pltype = {4 6 3 3};
title("Inference on inflation trend");
_plegctl={2 3.8 5.7 5};
_plegstr="Actual Inflation\000Trend Inflation\00090\% Confidence Bands";
xy(qtr[start:t],inf[start:t]~b_tt[.,1]~confint);

_pltype = {6 3};
_plegctl={2 4 6 5};
 _plegstr="Inflation Cycle\000Inflation Gap";
title("inflation cycle and inflation gap");
XY(qtr[start:t], xfor~inf[start:t]-b_tt[.,1]);
 
output file=UC1.out reset;

"== FINAL OUTPUT ===============================";
"Likelihood value is ";; -fout;
"Estimated parameters (s.e.) tstat are:"; 
prm_fnl~sd_fnl~prm_fnl./sd_fnl;
"===============================================";
output off;
end;

@ end OF MAIN PROGRAM @
@========================================================================@
 
PROC LIK_FCN(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
xfor0,xfor1,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2, 
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22;  

EXTERNAL PROC TRANS, V_PROB;
 
           PRMTR=TRANS(PRMTR1);

           P00=PRMTR[1,1];    @Pr[St=0/St-1=0]@
           P11=1;    @Pr[St=2/St-1=2]@
        
           siguGM0 = prmtr[2,1];
           siguGM1 = prmtr[3,1];
           sigv= prmtr[4,1]; 
           sigep = prmtr[5:6,1];
           rho = prmtr[7,1];
           phi1 = prmtr[8,1];
           phi2 = prmtr[9,1];
           psi=prmtr[10:11,1];    
           k=prmtr[12:13,1];

           A=1~0;
            
           H = 1~1~0~0|0~0~1~0;
                         
             
           B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
           P_LL0 = eye(4)*10000;          @uncertainty associated with initial guess@
           B_LL1=inf[t1]|0|xx[1]|0; 
           P_LL1 = eye(4)*10000;       
               
           subF = phi1~phi2| 1~0; 
           subQ= siguGM0^2~0|0~0;
           vecp = inv(eye(4)-subF.*.subf)*vec(subq);
           vv = reshape(vecp',2,2);
           p_ll0[3:4,3:4] = vv;
           p_ll1[3:4,3:4] = vv;
               
           p_ll0[2,2] = sigep[1,1]^2/(1-psi[1]^2);
           p_ll1[2,2] = sigep[2,1]^2/(1-psi[2]^2);
           
           PROB__0 = 1;
           PROB__1 = 0;
            

/*  START ITERATION */
/*===================================================================== */

       LIKV=0.0;
       J_ITER = 1;
       do until J_ITER>T;
    
            yt=inf[j_iter]|xx[j_iter];
     
            dmt = dmat[j_iter];
            sigu=siguGM0*dmt + siguGM1*(1-dmt);
 
            Q0=    sigv^2~0~0~0|
                   0~sigep[1]^2~rho*sigu*sigep[1,1]~0|
                   0~rho*sigu*sigep[1,1]~sigu^2~0|
                   0~0~0~0;
            Q1=    sigv^2~0~0~0|
                   0~sigep[2]^2~rho*sigu*sigep[2,1]~0|
                   0~rho*sigu*sigep[2,1]~sigu^2~0|
                   0~0~0~0;
            F0 =   1~0~0~0 |
                   0~psi[1]~0~0|
                   0~0~phi1~phi2|
                   0~0~1~0;
            F1 =   1~0~0~0 |
                   0~psi[2]~0~0|
                   0~0~phi1~phi2|
                   0~0~1~0;
            FF= phi1~phi2|1~0;
            xfor0 = k[1] *A*(inv(eye(2)-FF));           /*1x2 coef*/
            ome0 =0~0~xfor0| 0~0~0~0;
            xfor1 = k[2] *A*(inv(eye(2)-FF)); 
            ome1 =0~0~xfor1| 0~0~0~0;
           
        /* KALMAN FILTER */
       /*==============================================================*/
  
             B_TL00 = F0 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
             B_TL01 =F1 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@
             
             B_TL10 =F0 * B_LL1 ;  @WHEN S_{t-1}=1, S_{t}=0@
             B_TL11 = F1 * B_LL1 ;  @WHEN S_{t-1}=1, S_{t}=1@
         
             P_TL00 = F0* P_LL0 * F0' + Q0;
             P_TL01 = F1 * P_LL0 * F1' + Q1;
          
             P_TL10 = F0 * P_LL1 * F0' + Q0;
             P_TL11 = F1* P_LL1 * F1' +   Q1;
          
 
             F_CAST00= yt - H * B_TL00 - ome0* B_TL00;
             F_CAST01= yt- H * B_TL01 - ome1* B_TL01;
             
             F_CAST10= yt- H * B_TL10 - ome0*  B_TL10;
             F_CAST11= yt- H * B_TL11 - ome1* B_TL11 ;
            
 
             SS00= H * P_TL00 * H' ; @VARIANCE OF FORECAST ERROR@
             SS01= H*  P_TL01 * H' ;
            
             SS10= H * P_TL10 * H' ;
             SS11= H * P_TL11 * H' ;
        
 
             B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
             B_TT01 = B_TL01 + (P_TL01 * H' *inv( SS01)) * F_CAST01;
            
             B_TT10 = B_TL10 + (P_TL10 * H' *inv( SS10)) * F_CAST10;
             B_TT11 = B_TL11 + (P_TL11 * H' *inv( SS11) )* F_CAST11;
           

             P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
             P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
           
             P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
             P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
         
       /* HAMILTON FILTER */      /*==============================================================*/

             PR_VL00=V_PROB(F_CAST00,SS00)*p00 *PROB__0;    @Pr[St,Yt/Yt-1]@
             PR_VL01=V_PROB(F_CAST01,SS01)*(1-p00)*PROB__0;
           
             PR_VL10=0;
             PR_VL11=V_PROB(F_CAST11,SS11)*p11 *PROB__1;
         
                                          @CONDITIONAL DENSITY TIMES WEIGHT@

             PR_VAL=PR_VL00+PR_VL01+PR_VL10+PR_VL11;
                                  @WEIGHTED AVERAGE OF CONDITIONAL DENSITIES:
                                   f(y_t|Y_{t-1})@



             PRO_00=PR_VL00/PR_VAL;       @Pr[St,St-1/Yt]@
             PRO_01=PR_VL01/PR_VAL;
        
             PRO_10=PR_VL10/PR_VAL;
             PRO_11=PR_VL11/PR_VAL;
         
             PROB__0=PRO_00+PRO_10;       @Pr[St=0/Yt]@
             PROB__1=PRO_01+PRO_11;       @Pr[St=1/Yt]@
            
       /* COLLAPSING TERMS */
       /*==============================================================*/

              if(prob__0!=0);
              B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10 )/PROB__0;
              endif;
              if(prob__1!=0);
              B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11 )/PROB__1;
              endif;
               
              if(prob__0!=0);
              P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
                     PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)'))/PROB__0;
              endif;
              if(prob__1!=0);
              P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
                     PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)'))/PROB__1;
              endif;
              
              if J_ITER < START; goto skip; endif;
 
              LIKV = LIKV -LN(PR_VAL);
                skip:
       J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
/*===================*/


RETP(LIKV);
ENDP;
 
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
proc TRANS(prmtr1); @ constraining values of reg. coeff.@
local prm_new,d0,d1;
prm_new=prmtr1; 

    prm_new[1,1]=EXP(prmtr1[1,1]) ./ (1 + EXP(prmtr1[1,1])); 
    prm_new[10:11,1] =prmtr1[10:11,1]./(1+abs(prmtr1[10:11,1]));
    prm_new[7] =prmtr1[7,1]./(1+abs(prmtr1[7,1]));
    prm_new[2:6,1] = sqrt(prmtr1[2:6,1]^2);
    d0 =prmtr1[8,1]./(1+abs(prmtr1[8,1]));
    d1 =(1-abs(d0))*prmtr1[9,1]./(1+abs(prmtr1[9,1]))+ abs(d0)-d0^2;
    prm_new[8,1]=2*d0;
    prm_new[9,1]=-1*(d0^2 + d1) ;
  
retp(prm_new);
endp;


@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
PROC V_PROB(EV, HE);      @ CALCULATES    Pr[Yt/St,Yt-1]@

         local VAL;
         VAL=(1/SQRT(2*PI*det(HE)))*EXP(-0.5*EV'inv(HE)*EV);

RETP(VAL);
ENDP;
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@


proc (6) = FILTER(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
xfor0,xfor1,xfor2,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,confint,lowerCI,upperCI,stderr,stdprederr,xforecast,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2,B_M,PR_TTM,PR_TLM,twostep,
pr_tl00,pr_tl01, pr_tl02, pr_tl10,pr_tl11,pr_tl12,pr_tl20,pr_tl21,pr_tl22,
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22,
pi2hat_00,pi2hat_01,pi2hat_02,pi2hat_10,pi2hat_11,pi2hat_12,
pi2hat_20,pi2hat_21,pi2hat_22,B_T2L00,B_T2L01,B_T2L02,
B_T2L10,B_T2L11,B_T2L12,B_T2L20,B_T2L21,B_T2L22;  

EXTERNAL PROC TRANS, V_PROB;


           PRMTR=TRANS(PRMTR1);
 
            P00=PRMTR[1,1];    @Pr[St=0/St-1=0]@
           P11=1;    @Pr[St=2/St-1=2]@
        
           siguGM0 = prmtr[2,1];
           siguGM1 = prmtr[3,1];
           sigv= prmtr[4,1]; 
           sigep = prmtr[5:6,1];
           rho = prmtr[7,1];
           phi1 = prmtr[8,1];
           phi2 = prmtr[9,1];
           psi=prmtr[10:11,1];    
           k=prmtr[12:13,1];

           A=1~0;
            
           H = 1~1~0~0|0~0~1~0;
                         
             
           B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
           P_LL0 = eye(4)*10000;          @uncertainty associated with initial guess@
           B_LL1=inf[t1]|0|xx[1]|0; 
           P_LL1 = eye(4)*10000;       
               
           subF = phi1~phi2| 1~0; 
           subQ= siguGM0^2~0|0~0;
           vecp = inv(eye(4)-subF.*.subf)*vec(subq);
           vv = reshape(vecp',2,2);
           p_ll0[3:4,3:4] = vv;
           p_ll1[3:4,3:4] = vv;
               
           p_ll0[2,2] = sigep[1,1]^2/(1-psi[1]^2);
           p_ll1[2,2] = sigep[2,1]^2/(1-psi[2]^2);
           
           PROB__0 = 1;
           PROB__1 = 0;
            
           B_M=ZEROS(T,4);  @will store (1,1) element of B_TT@
           PR_TTM=ZEROS(T,2);  @will store Pr[S_t=0|Y_t]@
           PR_TLM=ZEROS(T,2);  @will store Pr[S_t=0|Y_{t-1}]@
        
           stdprederr=zeros(t,1);
           confint=zeros(t,2);  
           xforecast=zeros(t,1);
  
/*  START ITERATION */
/*===================================================================== */

       LIKV=0.0;
       J_ITER = 1;
       do until J_ITER>T;
             yt=inf[j_iter]|xx[j_iter];    
           
             pr_tl00 = PROB__0*p00;         
             pr_tl01 = PROB__0*(1-p00);
         
             pr_tl10 = 0;
             pr_tl11 = PROB__1*p11;
             PR_TLM[J_ITER,.]= (pr_tl00 +pr_tl10)
                                ~(pr_tl01 +pr_tl11 ) ; 
 
            dmt = dmat[j_iter];
            sigu=siguGM0*dmt + siguGM1*(1-dmt);
             
            Q0=    sigv^2~0~0~0|
                   0~sigep[1]^2~rho*sigu*sigep[1,1]~0|
                   0~rho*sigu*sigep[1,1]~sigu^2~0|
                   0~0~0~0;
            Q1=    sigv^2~0~0~0|
                   0~sigep[2]^2~rho*sigu*sigep[2,1]~0|
                   0~rho*sigu*sigep[2,1]~sigu^2~0|
                   0~0~0~0;
            F0 =   1~0~0~0 |
                   0~psi[1]~0~0|
                   0~0~phi1~phi2|
                   0~0~1~0;
            F1 =   1~0~0~0 |
                   0~psi[2]~0~0|
                   0~0~phi1~phi2|
                   0~0~1~0;
            FF= phi1~phi2|1~0;
            xfor0 = k[1] *A*(inv(eye(2)-FF));           /*1x2 coef*/
            ome0 =0~0~xfor0| 0~0~0~0;
            xfor1 = k[2] *A*(inv(eye(2)-FF)); 
            ome1 =0~0~xfor1| 0~0~0~0;
            
             B_TL00 = F0 * B_LL0 ;  @WHEN S_{t-1}=0, S_{t}=0@
             B_TL01 =F1 * B_LL0 ;  @WHEN S_{t-1}=0, S_{t}=1@
           
             B_TL10 =F0 * B_LL1 ;  @WHEN S_{t-1}=1, S_{t}=0@
             B_TL11 = F1 * B_LL1 ;  @WHEN S_{t-1}=1, S_{t}=1@
             

             P_TL00 = F0* P_LL0 * F0' + Q0;
             P_TL01 = F1 * P_LL0 * F1' + Q1;
              
             P_TL10 = F0 * P_LL1 * F0' + Q0;
             P_TL11 = F1* P_LL1 * F1' +   Q1;
             
             F_CAST00= yt - H * B_TL00 - ome0* B_TL00;
             F_CAST01= yt- H * B_TL01 - ome1* B_TL01;
            
             F_CAST10= yt- H * B_TL10 - ome0*  B_TL10;
             F_CAST11= yt- H * B_TL11 - ome1* B_TL11 ;
             

            xforecast[j_iter,.]=(xfor0*B_TL00[3:4,1])*pr_tl00 
                        +(xfor1*B_TL01[3:4,1])*pr_tl01 
                        +(xfor0*B_TL10[3:4,1])*pr_tl10 
                        +(xfor1*B_TL11[3:4,1])*pr_tl11 ;
  
             SS00= H * P_TL00 * H' ; @VARIANCE OF FORECAST ERROR@
             SS01= H*  P_TL01 * H' ;
           
             SS10= H * P_TL10 * H' ;
             SS11= H * P_TL11 * H' ;
        

             B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
             B_TT01 = B_TL01 + (P_TL01 * H' *inv( SS01)) * F_CAST01;
             
             B_TT10 = B_TL10 + (P_TL10 * H' *inv( SS10)) * F_CAST10;
             B_TT11 = B_TL11 + (P_TL11 * H' *inv( SS11) )* F_CAST11;
              
             P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
             P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
             
             P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
             P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
             
       /* HAMILTON FILTER */
       /*==============================================================*/


             PR_VL00=V_PROB(F_CAST00,SS00)*p00 *PROB__0;    @Pr[St,Yt/Yt-1]@
             PR_VL01=V_PROB(F_CAST01,SS01)*(1-p00)*PROB__0;
             
             PR_VL10=0;
             PR_VL11=V_PROB(F_CAST11,SS11)*p11 *PROB__1;
             
                                          @CONDITIONAL DENSITY TIMES WEIGHT@

             PR_VAL=PR_VL00+PR_VL01+PR_VL10+PR_VL11;
                                  @WEIGHTED AVERAGE OF CONDITIONAL DENSITIES:
                                   f(y_t|Y_{t-1})@

             PRO_00=PR_VL00/PR_VAL;       @Pr[St,St-1/Yt]@
             PRO_01=PR_VL01/PR_VAL;
             
             PRO_10=PR_VL10/PR_VAL;
             PRO_11=PR_VL11/PR_VAL;
            
            B_TT=B_TT00*PRO_00+B_TT01*PRO_01+B_TT10*PRO_10+B_TT11*PRO_11;
            B_M[J_ITER,.]=B_TT';
  
 
            lowerCI=(B_TT00[1,1]-1.645*sqrt(P_TT00[1,1]))*PR_tl00 
                    +(B_TT01[1,1]-1.645*sqrt(P_TT01[1,1]))*PR_tl01
                    +(B_TT10[1,1]-1.645*sqrt(P_TT10[1,1]))*PR_tl10
                    +(B_TT11[1,1]-1.645*sqrt(P_TT11[1,1]))*PR_tl11;
            upperCI =(B_TT00[1,1]+1.645*sqrt(P_TT00[1,1]))*PR_tl00 
                    +(B_TT01[1,1]+1.645*sqrt(P_TT01[1,1]))*PR_tl01
                    +(B_TT10[1,1]+1.645*sqrt(P_TT10[1,1]))*PR_tl10 
                    +(B_TT11[1,1]+1.645*sqrt(P_TT11[1,1]))*PR_tl11; 

            confint[j_iter,.] = lowerCI ~ upperCI;


            stderr= (f_cast00[1,1]/sqrt(ss00[1,1]))*pr_tl00 
                    +(f_cast01[1,1]/sqrt(ss01[1,1]))*pr_tl01
                    +(f_cast10[1,1]/sqrt(ss10[1,1]))*pr_tl10
                    +(f_cast11[1,1]/sqrt(ss11[1,1]))*pr_tl11;
  
            STDPREDERR[j_iter,.] = stderr;
            PROB__0=PRO_00+PRO_10;       @Pr[St=0/Yt]@
            PROB__1=PRO_01+PRO_11;       @Pr[St=1/Yt]@
          
            PR_TTM[J_ITER,.]=PROB__0~prob__1;

       /* COLLAPSING TERMS */
       /*==============================================================*/

              if(prob__0!=0);
              B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10 )/PROB__0;
              endif;
              if(prob__1!=0);
              B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11 )/PROB__1;
              endif;
             
              if(prob__0!=0);
              P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
                     PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)'))/PROB__0;
              endif;
              if(prob__1!=0);
              P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
                     PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)'))/PROB__1;
              endif;
               
              IF J_ITER < START; goto skip; endif;
              LIKV = LIKV -LN(PR_VAL);

                skip:
       J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
/*===================*/

retp(B_M[START:T,.],PR_TTM[START:T,.],PR_TLM[START:T,.],stdprederr[start:t,.],
         xforecast[start:t,.],confint[start:t,.]);
ENDP;


@==========================================================@

